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Predicting which crystalline modifications can be present in a chemical system requires the global 
exploration of its energy landscape. Due to the large computational effort involved, in the past this 
search for sufficiently stable minima has been performed employing a variety of empirical potentials 
and cost functions followed by a local optimization on the ab initio level. However, this entails 
the risk of overlooking important modifications that are not modeled accurately using empirical 
potentials. In order to overcome this critical limitation, we develop an approach to employ ab initio 
energy functions during the global optimization phase of the structure prediction. As an example, 
^ , we perform a global exploration of the landscape of LiF on the ah initio level and show that the 

O ' relevant crystalline modifications are found during the search. 
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I. INTRODUCTION 



A fundamental issue in solid state theory is the crystalline structure a given chemical system exhibits in the solid 
O . stateii^i^ii^i^. Why is a particular periodic atomic configuration adopted, which among several modifications is the 
' preferred one at a particular temperature and pressure, and which thermodynamically metastable but kinetically 
stable modifications are possible in the first place? Answering these questions requires the global exploration of the 
energy landscape of the chemical system^'*''^. Every metastable modification of a solid compound corresponds to 
a locally ergodic region on the energy landscape^, i.e. a set of atomic configurations which exhibits the property 
• that the equilibration time r^q of the chemical system within the region is much smaller than the observational time 
^ ' scale tabs, which in its turn is much smaller than the time scale Tesc on which the system is expected to leave this 
region, Teq <IC tobs ^ Tesc- In particular at low temperatures such locally ergodic regions constitute basins around 
one or several local minima of the energy landscape, and the kinetic stability of the corresponding compounds is 
controlled by the energetic and entropic barriers surrounding the region^i^. Thus, the first step in the prediction of 
the possible structures in a chemical system is the determination of the local minima on the energy landscape using a 
^ , global optimization algorithm. One should note that it is not sufficient to obtain only the global minimum: all local 
minima that are surrounded by sufficiently high energy barriers correspond to metastable modifications that may be 
^ of interest regarding their physical and/or chemical properties both in scientific and in technological applications^'^. 
^ Since the beginning of the 1990's, methods for theoretical structure determination and prediction employing global 

optimization techniques have been developed^ d°did^d'^d^d^ii6d''ii^d9i^°i^ii^^i^3 , using e.g. simulated annealing2ii2^, 
, genetic algorithm a^^i^^'^^ , or the threshold algorithm22i^. Recently, the combination of molecular dynamics with a 
C history dependent potential was suggested in the framework of the metadynamics approach^, in order to explore 

energy landscapes and phase transitions, e.gi^. The general idea, i.e. starting at a local minimum and exploring the 
^— I ' neighborhood on the landscape bears some resemblance to the lid or threshold algorith m??'?^ , but is different from 
T-H typical global optimization techniques such as simulated annealing, where long jumps on the energy landscape are 
r ■ allowed (sometimes called 'basin hopping'). The metadynamics employs molecular dynamics and focuses on a set 
' of a few relevant variables which are used to define excluded regions of the landscape and to describe the reaction 
J> , mechanism, while in the threshold algorithm the energy surface is stochastically sampled by a Monte-Carlo random 
walk below a sequence of energy thresholds, and the lid algorithm systematically explores the landscape below such 
. , energy lids by excluding all parts of phase space that have already been visited. 
' Since a typical set of global optimization runs involves millions or even billions of energy evaluations, a modular 
" ■ ' approach has become standard^^^, where a global search on an empirical energy / cost function landscape generates 
structure candidates, which are subsequently locally optimized on full quantum mechanical level using e.g. the Hartree- 
Fock approximation or density functional theoryiS.. Of course, this use of empirical potentials contains the risk that 
good candidates are overlooked because they do not correspond to a minimum (or only to a high-lying shallow one) 
on the empirical landscape, and there are many chemical systems where no straightforward empirical energy function 
based on simple or refined potentials or a crystal-chemically inspired cost function such as a bond-valence potential 
exist. But even for those systems such as ionic compounds where supposedly good model potentials have been 
constructed, the question to what degree the empirical energy landscape globally agrees with the ah initio energy 
landcape has been debated since the inception of work on structure prediction^. i'^^'^^ . Clearly, a careful comparison 
between these two energy landscapes for a particular system should yield much insight into the foundations of the 
current standard modular approach to structure prediction in solids. Only now the computers are reaching the speed 
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and ubiquity that will allow us to perform the global optimization on the ab initio level, as Oganov and co-workers2^ 
have shown using a hybrid genetic algorithm for this purpose. 

However, the step from model potentials to ab initio calculations is absolutely non-trivial, and requires careful 
adjustment of parameters to use only a miminum of CPU time and thus keep the calculations tractable. In this 
study, we investigate the energy landscape of the LiF-system using stochastic simulated annealing. This system was 
chosen as a test case, since we had studied the landscapes of the alkali halides in earlier workli using a Coulomb- plus 
Lennard-Jones-potential. Thus, the most important structure candidates in the system are known at the empirical 
potential level, and we can better judge both the success of the global exploration on the ab initio level and the degree 
of agreement between the empirical and ab initio energy landscapes than would have been possible when choosing a 
not-yet-investigated chemical system. 

There are thus two major goals of this article: Firstly, to show that the ab-initio exploration of energy landscapes 
with a Monte Carlo random walk based technique such as simulated annealing is feasible. Secondly, to investigate 
to what degree crucial features of the landscape such as the relevant minima on the level of the empirical potentials 
and on the ab initio level are the same. Finally, we note that the ab initio energy landscape is expected to be an 
appropriate choice for any system, whether ionic, covalent, or metallic. Thus, being able to globally explore such an 
ab initio energy landscape will open the path to structure prediction in systems which can no longer be reasonably 
described with straightforward empirical potentials. 

II. METHODS 
A. General approach. 

Our general approach to the determination of structure candidates has been given in detail elsewhere^i2^. To 
summarize: First, the minima on the energy landscape are determined using simulated annealing, possibly combined 
with a stochastic quench, as a global optimization algorithm, where both atom positions and cell parameters are freely 
varied. Next, the corresponding configurations are analyzed regarding their symmetries using an algorithm to find 
symmetrica^ and the space group^* as implemented in the program KPLOT'^^, followed by a comparison using an 
algorithm to compare cells^S,^ in order to eliminate duplicate structures. Finally, the structure candidates are locally 
optimized on the ab initio level using both a heuristic algorithm'*^''*^ and the energy minimization routines included 
in the various ab initio codes. We always employ several ab initio methods (Hartree-Fock and density functional 
theory), since this allows us to compare the ranking of the candidates by energy and thus to gain some estimate of 
their thermodynamic stability. This is crucial since, by definition, no comparison of the predicted structures with 
experimental data is possible and thus we cannot "tune" the parameters of the quantum mechanical methods to 
reproduce the experiment. Furthermore, comparing the outcomes of the local optimizations for different methods 
yields insights into the connectedness of the candidates via low-lying saddles on the energy landscape. Finally, if 
sufficient computational power is available, we can employ the lid or threshold algorith m^^i'^° , in order to quantitatively 
study the energetic and entropic barriers on the landscape, which control the kinetic stability of the metastable 
modifications of the solid compound. 

B. Ab initio calculations and global exploration: technical details. 

For the ab initio energy calculations we employ the program CRYSTAL2006^. A set of preliminary tests was 
performed to optimize the efficiency of our approach when applied to structure prediction on the ab initio level. The 
most important parameters tested were: the basis sets for the ab initio calculations; parameters such as integral 
thresholds for the ab initio calculations; the length of the simulated annealing run and of the subsequent quench run; 
the move classes involved. 

This preliminary step is actually crucial for the task of performing ab initio explorations of energy landscapes: the 
energy calculations are performed without the use of symmetry, since all possible structures must be accessible during 
the random walk. However, with the default parameters in CRYSTAL2006, a single Hartree-Fock calculation for an 
eight-atom simulation cell whose side length equals to the experimental lattice constant without symmetry (space 
group PI) takes ~ 13 minutes. A typical run may consist of 100000 and more simulated annealing steps, and thus 
the total CPU time for a single run would be on the scale of 10^ minutes, i.e. roughly 2 years. For the exploration of 
a landscape, dozens or even hundreds of such runs are necessary. In addition, it is necessary to achieve convergence 
of the self-consistency cycles in calculations which start from a random geometry. This makes it obvious that a very 
careful calibration of all parameters is necessary. One has to exploit the fact that only a rough knowledge about the 
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possible local mimima is required in the first stage of the global optimization. The final local optimization, based on 
an optimization via analytical gradients, can be subsequently done with good parameters. 

The initial tests resulted in the following choices: The basis sets from reference^ were selected, with a slightly 
modified fluorine basis set: slightly tighter sp functions were chosen for the two outermost exponents (0.45 instead 
of 0.437, 0.2 instead of 0.147) to enhance the numerical stability and the speed of the calculations. As the global 
optimization only has to provide a rough information about the energy landscape, it is not necessary to converge 
the solution very accurately. The threshold for the convergence of the self-consistent field (SCF) cycle was therefore 
reduced from 10~^ to 10~^ Eh- Similarly, the thresholds for neglecting integrals were reduced from the default values 
of 10-*', 10-^, 10-^, 10"^ 10^12 to 10-4^ 10-4^ 10-4^ 10-4^ 10-8. For the fc-point sampling, a shrinking factor of 2 
was used. The calculations were performed at the Hartree-Fock level. Note that since only the approximate positions 
of the basins need to be found during the initial global exploration, the level of theory does not play a crucial role. 

The initial cell was cubic with a cell parameter of 7.07 A. 4 lithium and 4 fluorine atoms were randomly placed 
in this cell. No symmetry was used, i.e. the simulated annealing and quenching was done in PI. The probabilities 
of the individual moves were as follows: moving individual atoms (70 %; maximal step size was 5 pm), exchanging 
atoms (10 %), moving atoms with (10%) and without (5%) simultaneous change of the unit cell, and the change of 
the origin (5 %; this move is important when the cell is subsequently truncated due to a change of the cell vectors). 
For all the moves which change the cell parameters, the probability of a suggested move to shorten the cell was set to 
70%, in order to speed up the shrinking of the cell. To avoid atoms coming too close, a minimum distance of the sum 
of the radii, multiplied by 0.8, was required. The radii were determined by using the MuUiken charges and linearly 
interpolating between the tabulated radii^ for the neutral atoms and the ions. With this choice of parameters, it 
turned out that only a very short simulated annealing run was required (5000 steps, with the starting and final 
temperatures corresponding to 1 eV and 0.9 eV, respectively), followed by a quench of ^ 5000 steps. The reason why 
the simulated annealing part could be kept so short is probably due to the fact that there are only two atom types 
and a very simple bonding mechanism involved. Clearly, when testing this approach with more demanding systems, 
exhibiting covalent bonds or a larger number of atom types, considerably longer runs are to be expected. 

We generated different initial atom positions by using different start values for the random number generator. 
About 70 runs were performed, of which about half converged to reasonable structure candidate, whereas the other 
half remained in energetically unfavorable situations such as very low densities or two-dimensional structures. This 
could certainly have been improved by performing longer simulated annealing runs, but only at a much higher 
computational cost. 

After the quench was finished, the space group of the configuration was analyzed with the program KPLOT— , 
using algorithms to identify the symmetry^^ and to find the space group A subsequent local optimization was 
performed with the CRYSTAL code, using analytical gradients for the nuclear positions^^'^^ and the unit cell^i^ 
and the full geometry optimization as implemented in the present release^^i^. As this local optimization is not too 
demanding in terms of CPU time and as a high accuracy is desirable, the integral thresholds and the threshold for 
SCF convergence were set to the default values, and a shrinking factor of 4 was used for the k net. Also, the original 
basis set as in^ was used. The optimization was performed both at the Hartree-Fock level and at the level of the 
local density approximation (LDA). The fully optimized structures were again analyzed with KPLOT. 

The computational effort was typically a few days for the simulated annealing and subsequent quench runs, and a 
few minutes up to one hour for the local optimization, on a single CPU of a standard PC. 

III. RESULTS 

The results of these optimizations are displayed in Tables 1 and 2. Eight promising low-energy structure candidates 
were found (shown in Figs. [1]- [31 generated with XCrysDen^): the rock salt structure as observed experimentally, 
the zincblende structure (sphalerite), the wurtzite structure, the so-called 5-5 structure^' (an ionic analogue to the 
Bk structure of hexagonal BN), the NiAs structure and three structures with space group 62, 7 and 36, denoted 
LiF(I), LiF(II) and LiF(III), respectively (see Table 1 regarding the fraction of runs that resulted in the various 
structure types). LiF(I) and LiF(II) consist of nets of LiF4-tetrahedra, with the first one containing narrow channels 
and the second one resembling a twisted sphalerite or wurtzite structure. Finally, LiF(III) consists of a network 
of LiFs square-pyramids. All these structure candidates had also been observed in earlier global searches in alkali 
halide systemSi^ using emipirical potentials consisting of a Coulomb-term and a van-der-Waals-term, and one should 
note that LiF(I), LiF(II) and LiF(III) are quite typical representatives of the higher- lying local minima. This result 
clearly demonstrates, that the global exploration of the energy landscape on the ab initio level is feasible and provides 
reasonable structures. 
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IV. DISCUSSION AND CONCLUSION 

Concerning the accuracy of the ab initio calculations, we note that at the Hartree-Fock level, the wurtzite structure 
exhibits the lowest energy, in contrast to the experimental observation that LiF is found in the rock salt type. The 
failure of the Hartree-Fock approach to account for the proper ground state may be attributed to the neglect of the 
van der Waals interaction which is important for alkali halides as was shown in^i^. Although the van der Waals 
interaction is not considered in the LDA either, the LDA performs better and predicts the rock salt structure as the 
ground state, perhaps due to LDA's inherent tendency to over-bind and favor higher coordinations. These observations 
are in good agreement with ab initio calculations by Cancarevic et al^i^ using Hartree-Fock and density functional 
theory, which showed that for LiF the functionals LDA found the rocksalt type as the minimum structure, while 
Hartree-Fock and the hybrid functional B3LYP found the wurtzite type, respectively. Regarding the thermodynamic 
stability of the various modifications as function of pressure, there is no change compared to the results of Cancarevic 
et al^i^, which were based on the global exploration of the empirical potential based enthalpy landscapes of LiF 
for ten different pressures ranging from -16 GPa to -1-160 GPa. Here, we only give the transition pressures based 
on the LDA functional: /3-BeO wurtzite at ~ -8.5 GPa, wurtzite 5-5 at ~ -5.5 GPa, and 5-5 — > rock salt 
at ~ -5.0 GPa. Possible high-pressure phases might be the NiAs-type or the CsCl-type according to some density 
functionals; however, even if these modifications were to become thermodynamically stable, in both instances the 
transition pressures would be very high (> 100 GPa) beyond the range of the validity of straightforward ab initio 
calculations. 

The comparison between the outcome of the global landscape explorations on the ab initio level and on the level of 
empirical potentials shows that the present investigation determines essentially all the minima with the lowest energies 
that had been found during the earlier work^^. Due to the relatively short simulated annealing runs, at the end of the 
global exploration phase many of the candidates had ended up in side-minima corresponding to distorted versions of 
the structures belonging to the main minimum of the various basins typically exhibiting lower symmetries. However, 
the subsequent local minimization using the standard high accuracy for the ab initio energy calculations resulted 
in the system reaching the main minima of the various basins. We note that in several instances the Hartree-Fock 
and the DFT (LDA) calculations reached different structure types, wurtzite and 5-5, respectively. This supports our 
earlier observations that the wurtzite and the 5-5-structure are close neighbors on the energy landscapci^ Similarly, 
we observed that the local mimization on DFT (LDA) level starting from four slightly distorted versions of the LiF (I) 
structure type, resulted in one case in the rock salt and in the other three cases in the 5-5-structure, respectively.— 
Again, this confirms our earlier results that the 5-5-type is located on the landscape close to the rock salt structure, 
possibly constituting a transitional modification on the route from the wurtzite to the rock salt type (c.f. Fig. [T|).— 
Finally, we note that in contrast to the study using the empirical energy landscape the structures exhibiting four- 
and five-fold coordination were found quite often. This is clearly a reflection of the fact that their energies on the ab 
initio level are very similar to the ones of the structure types with six- fold coordination (rock salt, NiAs) which are 
more strongly preferred when using the empirical potential. 

To summarize, the comparison between the empirical and ab initio energy landscapes shows that at least for the 
case of a simple ionic system such as LiF a) the minima representing the most relevant modifications and similarly 
most of the other chemically interesting structure types are present on both energy landscapes, b) however, their 
ranking in energy depends on the type of energy calculation and similarly the likelihood of observing the minima can 
also be quite different, c) structures which are closely related on the empirical landscape, i.e. separated by relatively 
small barriers, are also close neighbors on the ab initio energy landscape, and d) even classes of structures that are 
rather unusual such as those containing channels or square pyramids are found on both landscapes. We can thus 
conclude that one of the fundamental assumptions behind the standard approach of structure prediction of solids, i.e. 
that the empirical potentials can be employed in global landscape explorations to identify the relevant modifications 
of a solid compound is valid at least for ionic systems where reasonable suitable potentials are available. 

Of course, the number of simulated annealing runs possible when using the Hartree-Fock energy function is still 
much smaller than the number of runs with an empirical energy function. As a consequence, in particular the many 
high- lying and/or shallow minima associated with structures containing channels or belonging to transitions between 
two large basins, were detected relatively rarely. Nevertheless, the fact that all relevant modifications have been 
observed in the present study shows that the global exploration of the energy landscape of solids on the ab initio level 
using standard simulated annealing as the global optimization tool has finally become feasible. This will allow us to 
predict the possible modifications of crystalline compounds also in those chemical systems, where no simple empirical 
potentials can be constructed, thus overcoming one of the major hurdles facing crystal structure prediction in general 
chemical systems. 
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TABLE I: Structures obtained corresponding to local minima, the total energy (without a correction for the basis set superpo- 
sition error) for four formula units, in hartree units {Eh), and the number of times the structures were found. 

structure type space group energy (Eh) rmmber of times found 







HF 


LDA 


HF 


LDA 


rock salt 


225 


-428.2210 


-427.0665 


9 


11 


zincblende 


216 


-428.2211 


-427.0445 


10 


9 


5-5 


194 


-428.2222 


-427.0534 


4 


12 


wurtzite 


186 


-428.2250 


-427.0484 


8 


4 


NiAs 


194 


-428.2051 


-427.0515 


1 


1 


LiF(I) 


62 


-428.2162 




4 





LiF(II) 


7 


-428.2089 


-427.0374 


1 


1 


LiF(III) 


36 


-428.2054 




1 
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TABLE II; Optimized geometry of the various structures obtained. 
"Cell parameter from experiment (rock salt-type structure): a=4.027 A— 



structure type lattice constant, angle and Wyckoff positions 

(space group) HF LDA 

rock salt" (225) a=4.01 A a=3.94 A 

Li (0,0,0) (0,0,0) 

F (1/2,0,0) (1/2,0,0) 

zincblende (216) a=4.31 A a=4.23 A 

Li (0,0,0) (0,0,0) 

F (1/4, 1/4, 1/4) (1/4, 1/4, 1/4) 

5-5 (194) a=3.28 A; c=4.05 A a=3.24 A; c=3.94 A 

Li (1/3,2/3,1/4) (1/3,2/3,1/4) 

F (2/3,1/3,1/4) (2/3,1/3,1/4) 

wurtzite (186) a=3.09 A; c=4.86 A a=3.07 A; c=4.64 A 

Li (1/3,2/3,0) (1/3,2/3,0) 

F (1/3,2/3,0.386) (1/3,2/3,0.399) 

NiAs (194) a=2.79 A; c=4.80 A a=2.73 A; c=4.74 A 

Li (0,0,0) (0,0,0) 

F (1/3,2/3,1/4) (1/3,2/3,1/4) 

LiF(I) (62) a=5.60 A, 6=3.14 A, c=5.17 A 

Li (0.831, 3/4, 0.587) 
F (0.677, 1/4, 0.895) 

LiF(II) (7) a=2.92 A, &=5.62 A, c=5.34 A, (3 =115.9° a=2.84 A, 6=5.50 A, c=5.25 A; /3=115.5° 

Li (0, 0.407,0) (0, 0.405,0) 

Li (0.648, 0.121, 0.428) (0.638, 0.122, 0.430) 

F (0.955, 0.880, 0.676) (0.941, 0.878, 0.678) 

F (0.613, 0.403, 0.616) (0.599, 0.401, 0.617) 



LiF(III) (36) 



a=2.72 A, 6=5.30 A, c=4.81 A 
Li (0, 0.083, 0) 
F (0, 0.617,0.714) 
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FIG. 1: Ball-and-stick models of the structurally related rock salt (left), 5-5 (middle) and wurtzite (right) structures, viewed 
along the c-axis. Large and small spheres correspond to Li- and F-atoms, respectively. The wurtzite-type transforms into the 5- 
5-type by slightly displacing the Li- and F-atoms along the c-axis such that the local coordination changes from LiF4-tetrahedra 
to LiFs-trigonal bipyramids. Similarly, one obtains the rock salt-type from the 5-5-type by compressing the 5-5-structure along 
the a-axis such that the LiFs-trigonal bipyramids become LiFs-'square bipyramids', i.e. LiFs-octahedra. 




FIG. 2: Ball-and-stick models of LiF(I) (left), LiF(II) (middle) and LiF(III) (right), respectively. For notation, c.f. Fig. [T] 
These networks of LiF4-tetrahedra and LiFs-square pyramids are characteristic of metastable higher-lying local minima on the 
energy landscape of the alkali halides.-^'* 





FIG. 3: Ball-and-stick models of the sphalerite (left) and NiAs (right) structure, respectively. For notation, c.f. Fig. [T] 



